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ABSTRACT 

Recent  Experiments  on  the  TMX-U  device  at  LLNL  have  indicated  the 
possibility  of  a  low  frequency  (  at  «  )  drift  wave  being  driven  unstable  by 

the  injection  of  neutral  beams  in  the  thermal  barrier  region.  Nonlocal  linear 
theory  is  performed  for  this  instability  using  slab  geometry  and  the  electrostatic 
approximation.  Particle  simulations  are  used  to  verify  the  linear  theory  and  to 
recover  relevant  nonlinear  phenomena.  Important  nonlinear  effects  include:  (1) 
ExB  trapping  (2)  relaxation  of  the  beam  density  profiles  ,  and  (3)  weakly  non¬ 
linear  electrons  (  eS <t>/Te  =  0.1  ).  These  nonlocal  simulations  are  compared 
with  previous  simulations  of  Ref.  1.  which  used  a  local  approximation  for  the 
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1.  Introduction 

Drift  waves  have  long  been  an  area  of  major  activity  in  the  controlled  fusion  community. 
In  particular  particle  simulation  techniques  have  been  used  many  times  to  examine  some  drift 
waves  in  great  detail2-7.  In  this  paper  a  relatively  new  type  of  drift  wave  instability  is  exam¬ 
ined.  This  paper  is  motivated  by  some  recent  experiments  on  TMX-U  that  have  indicated  that 
electrostatic  drift  waves  may  be  driven  unstable  by  injection  of  neutral  beams  in  the  thermal 
barrier  region1.  The  linear  theory  and  the  particle  simulations  in  Ref.  1  were  performed  using 
the  local  approximation.  The  purpose  of  this  paper  is  to  relax  the  local  approximation  in  both 
the  theoretical  treatment  and  in  the  particle  simulations. 

Section  2  describes  the  model  equilibrium  which  uses  a  slab  model  with  a  uniform  mag¬ 
netic  field  together  with  the  electrostatic  approximation.  Section  3  describes  the  linear  nonlocal 
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analysis,  the  simulation  techniques,  and  the  linear  simulation  behavior.  Section  4  describes 
the  dominant  nonlinear  behavior  from  the  particle  simulations  and  includes  a  discussion  of  the 
dominant  effect,  which  is  ExB  trapping8.  Section  5  contains  a  brief  summary  of  the  major 
results. 

2.  Mode]  Equilibria 

The  geometry  for  the  model  is  given  in  Fig.  1.  The  equilibrium  density  gradient  is  per¬ 
pendicular  to  the  magnetic  field.  There  may  be  an  equilibrium  electric  field  in  the  x  direction. 
The  low  beta  approximation  is  made,  which  makes  it  possible  to  take  Bq  as  independent  of  the 
x  coordinate.  The  beam  velocity  parallel  to  the  magnetic  field  Ub  is  taken  to  be  much  less  than 
the  Alfven  speed  VA .  Under  these  conditions  for  the  beam  velocity  and  the  Alfven  speed  the 
electrostatic  approximation  may  be  employed.  Finally,  quasineutrality,  which  as  shown  in  Ref.  1 
is  valid  in  the  limit  wj,  »  fl2,  is  used  to  close  the  system  of  equations.  As  in  Ref.  1,  finite 
gyroradius  effects  are  ignored 

The  constants  of  the  motion  for  particles  are  the  guiding  center  position, 

X  -  x  +  vy/Clc-,  (1) 

the  particle  energy, 

e  =  (1/2  )mv2  +  q<t>(x)\  (2) 

and  the  component  of  the  velocity  parallel  to  the  magnetic  field, 

vz  .  (3) 

Equilibria  may  then  be  specified  using  the  constants  of  the  motion  to  form  model  distribution 

functions  for  the  ions  and  for  the  electrons. 

The  most  basic  distribution  function  is  used  for  the  electrons;  that  is  a  Maxwellian  distri¬ 
bution  in  particle  energy  and  a  multiplicative  factor  depending  upon  the  guiding  center  position. 
The  model  distribution  function  is  therefore  given  by 

Fe(Ar,€,v.)  -  ge(X)  exp|(-(l/2)mv2  -I-  e<t>(x))/Te J  .  (4) 

For  the  electrons  finite  gyroradius  effects  should  be  completely  negligible.  Therefore  the  equili¬ 
brium  spatial  distribution  function  for  the  electrons  is  given  by 
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ne Or)  -  A,  (x)exp|e<Mx)/7f  J  (5) 

where  Nf(x)  may  be  arbitrarily  specified. 

The  ions  consist  of  two  populations.  One  population  is  treated  as  a  cold  background  The 
second  population  consists  of  ions  counter  streaming  along  the  magnetic  field.  The  model  dis¬ 
tribution  function  for  either  population  may  be  given  by 


Fi(X,t,v.)  «*  g,(A')exp(-  {«  -  1/2 mvz2  -  e0  } /7~j)  (6) 

r 

X  exp [ — .5 (yy — v0)2/ v,J]+exp [— .5(v.+v0)2/v„2] 
where  v„  is  the  thermal  velocity  parallel  to  the  magnetic  field  and  7,  is  a  parameter  related  to 
finite  ion  gyroradius.  For  the  purposes  of  this  paper  ion  gyroradius  effects  are  ignored  Finite 
ion  gyroradius  effects  are  removed  by  taking  the  limit  7]~  0.  At  that  limit  the  term 
exp(— {«  —  \/2mv?  —  e0)/7i)  reduces  to  a  delta  function  at  the  minimum  value  for  the  perpen¬ 
dicular  kinetic  energy  at  the  given  value  of  *  (  which  is  the  definition  of  €0  ).  This  minimum 
value  of  perpendicular  kinetic  energy  corresponds  to  the  ExB  drift  velocity  due  to  any  equili¬ 
brium  electric  field  in  the  x  direction 


Since  finite  gyroradius  effects  are  being  ignored,  the  guiding  center  variable  X  reduces  to 


X  -  x  +  Vy(x)/tlC,  (7) 

where  v^Or)  is  the  ExB  drift  velocity  and  is  a  function  of  only  the  x  coordinate.  Therefore  the 

spatial  distribution  function  for  the  ions  becomes 


n,  (x)  —  N,  (x)  (8) 

where  N,  (x)  is  a  completely  arbitrary  function  and  represents  the  total  ion  denstty  including  the 

background  ions  and  the  beam  ions. 

Quasineutrality  then  determines  the  equilibrium  through  the  relation 

Ne(x)  exp (e<t>(x)/Te)  -=  A, Or)  .  (9) 

In  addition,  we  may  require  that  the  equilibrium  is  field  free,  that  is,  the  equilibrium  electric 

field  is  zero  everywhere.  This  requirement  eliminates  the  free  energy  associated  with  the  shear 

in  the  ExB  velocity.  In  this  paper  linear  theory  will  be  performed  only  for  model  equilibria  that 

are  field  free.  However,  some  simulation  work  is  included  for  the  case  where  an  equilibrium 
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electric  field  is  present. 

3.  Linear  Theory  and  Linear  Simulation  Results 

Section  3A  contains  the  linear  theory  for  the  inhomogeneous  slab  model  of  the  previous 
section.  Particle  simulation  techniques  are  used  in  order  to  recover  linear  behavior  as  described 
in  Sec.  3B.  The  results  obtained  by  linear  theory  are  compared  with  the  results  obtained  from 
the  particle  simulations. 

3A.  Linear  Theory 

The  eigenmode  equation  for  our  problem  can  be  obtained  from  the  drift  kinetic  equation 
since  w  «  fic,  and  all  finite  gyroradius  effects  are  being  ignored.  The  drift  kinetic  equation  is 
then  linearized  and  Fourier  transformed  for  a  mode  with  finite  values  for  both  ky  and  k:.  The 
resulting  eigenmode  equation  is  then  easily  obtained  from  Eq.  (10-166)  of  Ref.  9  by  summing 
over  the  background  ion  component  and  the  ion  beam  component.  Resonant  electron  effects, 
which  are  small  corrections  for  this  instability1  ,  are  ignored  and  only  Debye  shielding  is 
retained  for  the  electrons.  The  resulting  eigenvalue  equation  is  given  by 


(10) 


where 


(11) 


(12) 


and  G(ky.,k:)  represents  the  dispersion  relation  in  the  local  approximation  with  kx  =  0  and  is 
given  by 


61n  (Nc) 


(13) 


Z'(f+)  4-  Z'(f  J 


5 


k2kle 

Here  cjp,  is  the  total  ion  plasma  frequency,  ,NC  is  the  cold  background  ion  density,  Nb  is  the 
total  ion  beam  density,  8  —  Nb/(2NC)  represents  the  density  of  one  ion  beam,  k 2  «=  k2  +  k2, 
kpf  is  the  electron  Debye  length,  Z  is  the  usual  plasma  dispersion  function,  and 
|±  —  (to  ±  k:  £/6)/(V2v„)  with  Ub  the  average  beam  velocity  along  the  magnetic  field.  The 
quantities  <up, ,  Nc,  Nb  koe  and  8  are  all  functions  of  the  x  coordinate. 

In  order  to  solve  Eq.  (10)  specific  choices  for  the  density  profiles  and  boundary  condition 
must  be  made.  The  intent  of  our  calculations  is  only  to  include  the  dominant  effects  of  the 
problem  in  the  most  straight  forward  manner  possible.  Therefore,  for  the  calculations  given  in 
this  paper,  we  choose  the  boundary  conditions  d<t>/dx  =  0  at  both  boundaries  The  ion  density 
profiles  are  given  the  form 

N,  Or)  =  nod  +  ^cos^x/L*))  (14) 

or 

N,(x)  =  rioO  +  B  sininx  /  Lx))  (15) 

and  may  vary  for  the  background  ions  and  the  beam  ions.  The  electron  density  is  then  deter¬ 
mined  from  quasineutrality;  that  is 

N,(x)  -  'ZN^Hx)  .  (16) 

5 

Equation  12  was  solved  by  integrating  the  equation  from  both  boundaries  and  matching 
the  results  in  the  center  of  the  system.  In  this  fashion  it  is  possible  to  recover  the  family  of 
complex  eigenmodes  {<f>ky  ^Gr)}  and  the  corresponding  eigenvalues  M.  For  those  cases  where 
the  local  approximation  is  valid  kyLn»  1  (  where  Ln  is  the  background  ion  density  scale 
length  )  the  local  approximation  gives  good  agreement  with  the  most  unstable  eigenvalue.  How¬ 
ever,  in  place  of  a  single  value  for  cu  a  set  of  eigenvalues  is  obtained.  This  fact  makes  it  more 
difficult  to  isolate  linear  behavior  in  the  simulation  code,  as  will  be  discussed  in  Section  3B 
For  example,  Fig  2.  shows  the  most  unstable  eigenmode  and  Fig  3  shows  the  second  most 
unstable  eigenmodes  for  a  test  case  with  the  parameters;  «=  0.03,  kvps  *=  0.74, 

k-L„  -=0.16,  Ub/c5  *  3.0,  \ulUb  *  0.0  ,  and  6 1 L  l2  =  0.1.  The  wave  vector  was  chosen  so  that 
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the  growth  rate  is  close  to  the  maximum  growth  rate  for  the  given  physical  parameters.  Here  c5 
is  the  sound  speed,  ps  L„  the  the  minimum  value  for  the  background  ion  scale 

length,  and  -  cj2L„  represents  the  maximum  drift  wave  frequency.  The  density  profile 
from  Eq.  (14)  is  used  with  Ab  —  —1.0  and  Ac  —  0.5.  The  growth  rates  are  very  close  to  each 
other  and  therefore  both  may  be  expected  to  be  present  is  simulations  with  random  excitation. 

The  particular  eigenmodes  shown  in  Fig.  2  and  Fig.  3  have  their  real  and  imaginary  com¬ 
ponents  almost  equal.  This  fact  means  that  the  eigenmode  has  approximately  elliptical  potential 
contours  in  the  plane  of  the  electric  field  with  one  axis  parallel  to  the  x  axis  .  This  type  of 
eigenmode  structure  is  retained  even  when  thermal  effects  are  added. 

Changing  the  beam  ion  density  profile  to  that  in  Eq.  (15)  with  the  value  Bb  =  1.0  yields 
eigenmodes  of  a  different  nature  as  seen  in  Fig.  4.  For  this  type  of  eigenmode  the  potential 
contours  may  be  considered  elliptical  about  the  potential  extrema  but  now  neither  axis  is  paral¬ 
lel  to  the  x  axis.  Thus  virtually  any  type  of  mode  structure  is  possible  for  this  problem. 

Keeping  ky  fixed  and  varying  Lx  causes  the  validity  condition  for  the  local  approximation 
to  be  invalid  and  it  is  expected  that  ymax/(D^ ax  would  decrease.  This  is  indeed  what  happens. 
However,  fif,/cumax  =  kyLn:  therefore  a  decrease  in  Lx  is  leads  to  a  larger  growth  rate.  In  this 
sense  finite  length  effects  do  not  serve  to  reduce  the  growth  rate. 

3B.  Simulation  Techniques  and  Linear  Behavior 

This  section  describes  the  simulation  techniques  and  the  linear  phenomena  recovered 
from  the  simulations.  A  quiet  start  technique  is  necessary  to  initialize  the  simulations  and 
achieve  good  agreement  with  the  theoretically  calculated  mode  structures  from  the  previous 
section. 

The  simulation  code  is  a  2  -  1/2  dimensional  electrostatic  particle  code  with  Boltzmann 
electrons.  The  field  solve  is  quasineutrality  as  in  Ref.  1.  The  geometry  is  that  of  Fig.  1  with  the 
simulation  plane  at  a  small  angle  to  the  x  —  y  plane  so  as  to  include  a  finite  value  for  the  paral¬ 
lel  wavevector  k..  The  code  employs  periodic  boundary  conditions  for  the  variation  perpendicu¬ 
lar  to  the  x  axis.  The  code  reflects  the  particles  at  x  —  0  and  x  —  Lx.  The  electric  field 
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component  Ex  is  zero  at  x  —0  and  x  —  Lx.  Full  ion  dynamics  are  retained.  This  has  the  disad¬ 
vantage  of  the  introduction  of  ion  cyclotron  waves.  However,  since  only  the  guiding  center 
motion  is  important  the  time  step  does  not  have  to  accurately  represent  the  high  frequency 
waves  present.  Only  a  single  beam  will  be  unstable  since  only  a  single  value  of  <o/k.  is  possible 
using  this  2  —  1/2  dimensional  model. 

In  order  to  recover  the  linear  behavior  with  the  nonlinear  particle  code  numerical  noise 
must  be  suppressed  as  much  as  possible  in  the  choice  of  initial  conditions.  Therefore,  for  the 
purposes  of  comparing  the  code  results  with  the  results  of  linea*-  theory  we  choose  a  case  where 
the  initial  thermal  spread  in  vz  is  zero.  In  addition,  the  particles  were  loaded  in  a  with  very  spe¬ 
cial  initial  conditions.  For  the  equilibrium  the  particles  were  loaded  uniformly  in  both  the  x 
direction  and  in  the  y  direction.  This  uniform  loading  guarantees  equivalent  particle  statistics  at 
each  location  in  the  simulation.  In  order  to  represent  the  density  gradient  it  is  necessary  to 
weight  the  particles  according  to  the  exact  value  of  the  density  at  the  position  they  are  loaded. 

These  techniques  are  ,  in  general,  not  enough  As  was  mentioned  in  the  previous  section 
the  generic  equilibrium  has  many  closely-spaced  eigenvalues.  Even  with  the  low  noise  level  in 
the  equilibrium  the  most  unstable  mode  will  not  be  able  to  become  dominant  before  the  simu¬ 
lation  becomes  nonlinear.  Therefore  special  initial  perturbations  are  required. 

One  approach  is  to  integrate  the  particles  in  the  fields  of  the  most  unstable  waves.  This 
may  be  done  either  numerically  or  analytically  When  this  was  done  either  way  the  self  con¬ 
sistent  field  as  computed  from  the  particles  was  very  close  to  the  theoretical  field  This  does 
suggest  that  the  calculated  solutions  are  indeed  self-consistent  solutions. 

By  using  this  type  of  initial  condition  the  simulation  can  then  be  carried  out  in  a  self- 
consistent  manner.  Items  of  interest  include  the  real  frequency,  the  growth  rate,  and  the  mode 
structure  in  the  inhomogeneous  x  direction.  The  physical  parameters  are  those  of  Fig.  2  and 
Fig.  3.  The  simulation  parameters  are  ftc,  At  -  0.5  -  1.0,  the  number  of  grids  in  the  x  direc¬ 
tion  is  Nx  -  65,  the  number  of  grids  in  the  periodic  direction  is  32,  and  roughly  six  thousand 
particles  for  the  background  ions  and  sixteen  thousand  particles  for  each  of  the  beams.  Larger 
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time  steps  can  be  taken  by  use  of  a  first  order  particle  integration  scheme  to  remove  the  time 
constraint  caused  by  the  high  frequency  waves10.  In  addition,  only  one  mode  was  kept  in  the 
periodic  direction  and  only  eight  modes  were  kept  in  the  inhomogeneous  x  direction.  The 
equilibrium  has  no  initial  electric  field. 

Figure  5  shows  the  the  potential  as  a  function  of  time  for  a  location  near  the  center  of 
the  mode  structure.  Notice  that  due  to  the  approximate  normal  mode  initialization  the  perturba¬ 
tion  begins  growing  with  a  definite  frequency  almost  immediately;  i.  e.  there  is  only  a  small 
component  present  that  does  not  comprise  part  of  the  normal  mode.  The  real  frequency  and 
growth  rate  are  within  about  5%  of  the  numerically  calculated  values.  Figure  6  shows  the 
recovered  eigenmode  from  the  simulation  during  the  linear  stage  using  the  Zed  analysis  pack¬ 
age11. 

The  results  given  in  Fig.  5  and  Fig.  6  are  good.  However,  the  agreement  is  certainly  not 
exact.  Sources  of  discrepancy  include  (1)  retaining  only  the  low  frequency  response  (  i.e. 
ignoring  terms  of  order  <o/Clci  ~  0.025  ),  (2)  not  including  all  of  the  Fourier  components 
necessary  for  complete  representation  of  the  mode  structure  (3)  the  excitation  of  cyclotron 
modes  and  (4)  the  excitation  of  secondary  normal  modes  due  to  the  sensitivity  of  the  initial 
excitation.  These  last  three  effects  will  now  be  discused  individually. 

Fourier  decomposition  of  the  calculated  mode  structures  allows  one  to  determine  how 
many  modes  to  keep  in  the  simulations.  Unfortunately  the  mode  structures  typically  have  a  tail 
of  a  few  percent  at  wavelengths  of  10  -  15  times  the  primary  wavelength. 

These  shorter  wavelengths  give  high  cyclotron  noise  unless  large  numbers  of  particles  are 
being  used  to  represent  the  ion  beams.  This  noise  is  due  primarily  from  the  extra  degree  of 
freedom  in  ion  phase  space  associated  with  the  gyromotion.  Thus  the  particle  motion  need  not  ( 
and  in  general  does  not )  contain  just  the  low  frequency  drift  motion.  Since  ui  «  flCI  the  high 
frequency  ion  cyclotron  waves  may  randomize  portions  of  the  particle  phase  space  during  an  e- 
folding  time  of  the  instability. 
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The  initial  excitation  requires  accurate  knowledge  of  the  real  frequency  and  the  growth 
rate  Deviation  from  the  true  values  yields  phase  errors  between  the  perturbed  positions  and 
velocities.  This  then  excites  other  normal  modes  with  a  finite  amplitude.  Some  of  these  other 
normal  modes  have  comparable  real  frequencies  and  growth  rates,  which  makes  interpretation 
of  simulation  results  more  difficult. 

4.  Nonlinear  Behavior 

This  section  describes  the  dominant  nonlinear  effects  as  observed  in  the  simulations, 
These  effects  include  (1)  ExB  trapping  of  resonant  beam  ions,  (2)  modification  of  the  density 
profiles  of  the  different  species,  and  (3)  weakly  nonlinear  electrons 

4 A  Nonlinear  Effects 


Some  of  these  effects  were  discussed  in  an  earlier  paper1  and  the  results  are  summarized 
in  Table  1.  These  results  were  obtained  by  assuming  a  local  approximation  for  the  fields,  that  is, 
that  the  perturbation  electric  field  is  assumed  to  be  zero  in  the  x  direction. 

We  wish  to  review  some  of  the  analysis  for  ihe  case  of  a  single  wave  <t>*~<t>0cos(kyy  +  k.z) 
in  the  local  approximation;  i.e.  there  is  no  x  dependence  in  the  problem.  The  approximate 
scaling  law  for  saturation  by  particle  trapping  in  the  local  approximation  is  given  by 


g&d>  __  y2  -> 

T,  k;W 


(17) 


which  is  obtained  by  equating  the  linear  growth  rate  of  the  instability  y  to  the  trapping  fre¬ 
quency  of  the  particle  in  the  wave  field  given  by  k.-JeSQ/m  .  (  The  factor  3  represents  an  "e- 
folding"  from  the  onset  of  nonlinear  effects  to  saturation  and  is  very  qualitative  in  nature  ) 
Associated  with  the  trapping  in  vz  is  displacement  in  the  x  direction  governed  by  the  relation 


v,  —  I  Oc,jr  —constant  .  (18) 

Equation  (18)  is  obtained  by  using  the  guiding  center  equations  of  motion.  Equation  (18) 
implies  that  trapped  particle  are  transported  in  x.  The  signs  of  the  various  quantities  in  Eq  (18) 
are  such  that  the  trapped  ions  (  particles  that  lose  energy  on  the  average  )  have  a  net  drift  up 
the  background  ion  density  gradient. 
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One  of  the  goals  of  this  paper  is  to  relax  the  local  approximation.  The  analysis  and  simu¬ 
lations  will  be  limited  to  the  case  of  a  single  value  for  (ky,k.).  As  mentioned  before,  only  one 
ion  beam  will  be  unstable  since  only  one  value  for  <o/k,  is  allowed  in  the  simulation  model.  An 
important  change  in  the  nonlocal  problem  over  the  local  problem  is  that  the  local  treatment 
assumes  that  there  is  no  perturbed  electric  field  parallel  to  the  density  gradient.  Particle  orbits 
are  constrained  so  that  their  ExB  drift  velocity  perpendicular  to  the  density  gradient  is  zero. 
Inclusion  of  the  correct  electric  field  means  that  particles  may  have  a  component  of  their  drift 
velocity  perpendicular  to  the  density  gradient.  A  second  important  feature  is  that  simulations 
using  the  local  approximation  employ  density  profiles  that  do  not  change  (  since  the  electric 
fields  are  independent  of  the  x  coordinate  ).  This  is  not  true  for  nonlocal  simulation  models. 

One  interesting  fact  is  that  the  constant  of  motion  given  in  Eq.  (18)  is  unchanged  even 
for  the  nonlocal  problem  as  long  as  there  is  only  a  single  value  for  the  wavevector  (ky,k.)\  that 
is,  the  potential  has  the  form 

<f>(xy,z,t)  *  <f>0(x,t)cos(kyy  +  kzz  -  ait)  (19) 

and  the  time  variation  of  the  potential  satisfies  d/dr  «  fif,  so  that  the  guiding  center  approxi¬ 
mation  may  be  made.  The  polarization  drift  velocity  is  very  small  compared  to  the  ExB  drift 
velocity  and  is  therefore  neglected  in  the  derivation  of  Eq.  (18). 

As  mentioned  above,  one  effect  missing  entirely  in  the  local  approximation  is  the  effect  of 
electric  field  parallel  to  the  density  gradient.  This  may  cause  ExB  trapping.  Idealizing  the 
problem  to  the  case  where  the  electric  field  is  entirely  perpendicular  to  the  magnetic  field  allows 
a  clear  intuitive  picture  of  ExB  trapping.  Consider  a  resonant  particle  moving  in  the  plane  per¬ 
pendicular  to  the  magnetic  field  and  retain  only  the  ExB  drift  velocity.  The  particle  attempts  to 
stay  on  equal  potential  contours  and  to  drift  around  them.  If  the  particle  may  drift  around  the 
potential  contours  and  change  its  position  significantly  (  comparable  to  the  scale  lengths  of  the 
problem  )  ,  then  linear  theory  is  no  longer  valid. 

An  estimate  for  when  ExB  trapping  is  important  can  be  obtained  from  equating  the  wave 
frequency  in  the  particle  frame  with  the  frequency  for  a  particle  to  drift  around  one  of  the 
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potential  peaks.  This  yields 


W  “  &>ExB 


(20) 


with 


wexb  =  (***>■)  ~  c?mci 


(21) 


where  w  is  the  wave  frequency  in  the  particle  frame.  The  condition  for  importance  of  this 
effect  is  then  given  by 


(22) 


The  analysis  is  actually  far  mo>-e  complicated  than  is  suggested  above  The  reason  is  that 
the  wave  has  a  nonzero  electric  field  component  in  all  three  important  directions  in  this  prob¬ 
lem.  The  effects  of  the  constant  of  the  motion  given  in  Eq.  (18)  must  be  combined  with  the 
motion  produced  by  Ex  in  determining  actual  particle  orbits.  Equation  (18)  implies  that  particles 
with  circular  type  motion  in  the  x  -  ^  plane  have  their  vr  changing  in  a  like  manner.  The  par¬ 
ticle  orbits  are  actually  very  problem  dependent  since  the  orbits  depend  on  the  phase  structure 
of  the  eigenmodes.  The  phase  structure  of  the  eigenmode  can  varv  greatly,  as  shown  in  Sec. 
2A. 

4B  Simulation  Results 

Nonlinear  phenomena  recovered  from  the  simulations  are  presented  in  this  section.  Con¬ 
nections  will  be  made  to  the  nonlinear  phenomena  discussed  in  the  previous  section;  (1)  ExB 
trapping,  (2)  the  final  appearance  of  the  density  profiles  and  (3)  saturation  levels.  An  exhaus¬ 
tive  parameter  search  is  entirely  out  of  the  question.  Parameters  that  will  be  varied  are  (1)  the 
thermal  spread  of  the  beams,  (2)  the  density  profile  of  the  ion  beams  and  (3)  whether  the 
equilibrium  has  a  zero  order  electric  field 

One  set  of  simulations  was  performed  in  order  to  investigate  the  effects  of  the  thermal 
spread  along  the  magnetic  field  on  the  saturation  of  the  instability.  The  physical  parameters  are 
those  corresponding  to  Fig.  2.  and  Fig.  3  except  that  Ab  -  0.0.  The  simulation  parameters  are 
A t  «=  0.5,  Nx  -  65  grids  in  the  x  direction  and  32  grids  in  the  periodic  direction  and  roughly 
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thirty  thousand  particles.  The  values  for  the  length  in  the  periodic  direction  was  chosen  so  that 
( ky,k. )  corresponded  to  the  values  for  the  most  unstable  root  in  the  local  approximation.  Eight 
modes  were  kept  in  the  inhomogeneous  direction.  The  thermal  spread  was  varied  from 
v,,/ Ub  «=  0.0  to  v„ / Ub  —  0.25.  These  were  the  same  simulation  parameters  used  for  previous 
simulations  using  the  local  approximation'.  Typically  many  orders  of  magnitude  in  growth  in 
the  mode  energy  was  observed  in  the  simulations  as  shown  in  Fig.  7. 

The  saturation  values  are  shown  in  Table  2  for  the  simulations.  As  can  be  seen  the  satura¬ 
tion  level  is  relatively  independent  of  the  beam  thermal  spread  parallel  to  the  magnetic  field 
unlike  the  simulation  in  Ref.  1.  This  result  is  very  suggestive  of  ExB  trapping  as  opposed  to 
trapping  in  the  direction  parallel  to  the  magnetic  field  since  the  saturation  lelel  does  not  depend 
much  on  the  growth  rate,  in  conflict  to  Eq.  (17).  Had  both  beams  been  unstable  (  which  would 
require  a  fully  three  dimensional  simulation  code  )  then  one  would  expect  that  the  peak  poten¬ 
tial  energies  would  be  on  the  order  of  20%  of  the  electron  thermal  energy. 

The  estimate  for  when  a  resonant  particle  would  become  nonlinearly  dominated  by  the 
motion  perpendicular  to  the  magnetic  field  is  approximately  when  the  ExB  trapping  frequency 
is  comparable  to  the  residual  wave  frequency  induced  by  parallel  trapping.  Expressing  this  in 
terms  of  an  approximation  yields 

kypskx\e<}>ITe^  =  Ar,|ed>/7'(,J1  2  (23) 

or  rewritten  as 


e<j>/Te  = 


k. 

2 

l 

A 

kyPs 

This  saturation  level  is  approximately  independent  of  the  growth  rate  for  this 
although  a  slight  dependence  still  exists  because  the  percent  of  resonant  particles 
in  determination  of  the  final  saturated  potential. 


(24) 

mechanism, 
plays  a  role 


More  evidence  for  the  saturation  mechanism  may  be  presented  by  a  careful  examination 
of  test  particle  orbits.  In  particular  particle  trajectories  in  the  plane  perpendicular  to  the  mag¬ 
netic  fields  are  informative.  A  test  particle  trajectory  for  one  simulation  is  shown  in  Fig.  8. 
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Notice  that  the  orbit  is  consistent  with  ExB  trapping;  that  is  particle  motion  in  both  directions 
perpendicular  to  the  magnetic  field  is  significant  Figure  9  shows  the  potential  as  a  function  of 
space  near  saturation  of  the  instability. 

Particle  trajectories  may  have  a  different  appearance  the  trajectory  shown  in  Fig.  8.  Figure 
10  gives  another  example  of  ExB  trapping  for  a  different  situation.  Notice  that  the  particle 
motion  in  both  directions  perpendicular  to  the  magnetic  field  is  significant  although  the  exact 
particle  motion  is  quite  different  from  that  in  Fig  9.  The  particle  motion  depends  upon  how 
close  that  particle  is  to  resonance  with  the  wave  and  upon  the  spatial  mode  structure  of  the 
wave.  Also,  even  with  a  single  value  for  \ky,kz),  several  eigenmodes  may  be  present.  Thus  the 
exact  motion  is  very  complicated  in  general. 

In  order  to  examine  the  constant  of  motion  given  by  Eq.  (19)  a  group  of  test  particles 
with  the  same  initial  values  of  x  and  vz  was  followed.  Constancy  of  the  entity  in  Eq.  (18) 
implies  that  the  particles  should  fall  on  a  line  in  x  —  vz  plane.  This  was  tried  and  proved  to  be  a 
useful  diagnostic  as  presented  in  Fig.  11.  Notice  that  particles  that  lost  energy  suffered  dis¬ 
placements  up  the  background  ion  density  gradient.  Since  more  ion  beam  panicles  lose  energy 
than  gain  energy  a  net  displacement  should  occur  in  the  ion  beam  density  profiles.  This  charac¬ 
teristic  is  seen  in  Fig.  12. 

The  displacement  of  the  ion  beam  need  not  be  small.  To  examine  this  a  very  unstable 
case  is  chosen  with  Ab  —  -1.0,  v„/Ub  *=  0.19  and  the  other  parameters  the  same  as  the  previ¬ 
ous  simulations.  This  type  of  ion  beam  density  profile  is  an  idealization  of  the  situation  that 
may  occur  when  the  injected  neutral  beam  does  not  penetrate  far  into  the  plasma.  Figure  13 
displays  the  changes  in  the  density  profiles  occurring  for  the  simulation.  The  situation  may  be 
clarified  further  by  examination  of  the  particle  distribution  functions  over  different  portions  of 
the  simulation  region  Those  regions  that  suffered  a  net  density  decrease  did  so  at  the  expense 
of  the  lower  energy  beam  ions.  Those  particles  were  transported  up  the  background  ion  density 
gradient.  These  considerations  are  seen  in  Fig.  14. 

Finally,  a  short  discussion  is  included  on  the  effects  caused  by  an  equilibrium  electric 
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field,  First  consider  the  case  of  a  uniform  electric  field  in  the  x  direction.  Formal  substitution 
into  the  local  or  nonlocal  dispersion  relation  shows  that  the  growth  rates  and  wavevectors  do 
not  change.  However,  the  real  frequencies  become 

t0±  =  ±A :zUb  —  ky  VE  (25) 

where  VE  —  Eq/B  is  the  zero  order  drift  velocity  in  the  y  direction.  For  typical  parameters  of 

this  problem  kyVE  is  not  negligible  to  k.  Ub .  For  the  case  where  the  electric  field  is  pointing 
toward  the  walls,  the  signs  of  the  quantities  in  Eq.  (25)  are  such  as  to  decrease  the  magnitude 
of  the  frequency  in  the  laboratory  frame.  In  the  extreme  limit,  this  effect  may  cause  the  per¬ 
pendicular  phase  velocity  to  be  in  the  ion  diamagnetic  drift  velocity  rather  than  in  the  electron 
diamagnetic  drift  direction. 

When  VE  is  spatially  dependent  the  situation  is  quite  different.  The  phase  velocity  of  the 
wave  in  the  frame  of  the  background  plasma  then  depends  upon  the  location  in  x.  Therefore, 
in  principle,  any  portion  of  the  distribution  function  may  be  in  resonance.  In  addition,  one 
might  expect  that  spatial  dependence  of  equilibrium  quantities  such  as  the  electron  tempera¬ 
ture  the  magnetic  field  and  the  zero  order  electric  field  would  cause  the  saturation  to  be  more 
mild.  Simulation  results  tend  to  support  this  idea.  An  example  is  shown  in  Fig.  15  which  shows 
the  change  in  the  linear  frequency.  Figure  16  shows  the  density  profiles  at  saturation. 

5.  Conclusions 

In  this  paper  both  linear  and  nonlinear  properties  of  an  ion-beam  driven  drift  instability 
are  described  using  an  inhomogeneous  slab  model.  The  agreement  between  linear  theory  and 
the  simulation  model  is  modestly  good.  The  linear  effects  of  some  of  the  parameters  are  dis¬ 
cussed.  In  addition,  the  particle  simulations  give  insight  into  the  nonlinear  saturation  of  the  ins¬ 
tability.  The  most  relevant  nonlinear  effects  for  present  experiments  are  (1)  ExB  trapping  of 
the  beam  ions,  (2)  cross  field  transport  of  the  trapped  beam  ions  and  (3)  weakly  nonlinear  elec¬ 
trons.  It  is  interesting  to  note  that  the  beam  ions  are,  on  average,  transported  up  the  back¬ 
ground  ion  density  gradient.  This  fact  suggests  that  the  instability  presented  here  may  actually 
serve  to  improve  confinement  of  the  high  energy  beam  ions. 
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Saturation  Characteristics 

====“==  m  = 

8 

small 

large 

weak 

ion  beam  trapping;  enhanced 

trapping  of  resonant  ions; 

cross-field  motion  of  trapped 

enhanced  cross-field  motion 

ions;  coherent  mode  coupling 

of  trapped  ions 

strong 

nonlinear  electrons;  non¬ 

trapping  of  resonant  ions; 

linear  motion  parallel  to  the 

enhanced  cross-field  motion 

density  gradient;  ion  beam 

of  trapped  ions;  weakly  non¬ 

trapping;  enhanced  cross  field 

linear  electrons 

motion  of  trapped  beam  ions 

TABLE  1.  Nonlinear  properties  as  a  function  of  parameter  space. 
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Saturation  Dependence  ot  thermal  Spread 

Vtf  /  ub 

- 4 — * - 

u.uu 
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cm 
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\cbipj Te)  sat 

“D.l  2 

0.12 

'  o.n — 

OTTO 

“0.082 

TABLE  2.  Variation  of  saturation  with  the  parallel  thermal  spread  of  the  ion  beams.  The 
very  weak  dependence  of  the  saturation  on  the  growth  rate  suggests  ExB  trapping. 
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FIG.  1.  Slab  geometry  for  the  instability.  The  magnetic  field  is  uniform  in  space. 
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complex  frequency  is  given  by  «,/«««  -  0.78  and  y/«w  0.  . 


4 


^ky,  k2M 


FIG.  3.  Numerically  obtained  solutions  for  the  second  most  unstable  solution  showing 
(a)  the  real  and  imaginary  parts  and  (b)  the  square  of  the  amplitude.  The  complex  fre¬ 
quency  is  given  by  (ur/wma*  “  0.74  and  y/fc>max  “  0.25. 
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FIG.  3b. 


FIG.  4.  Numerically  obtained  solutions  using  the  beam  ion  density  profile  from  Eq. 
(15)  showing  the  real  and  imaginary  parts  for  (a)  the  most  unstable  mode  and  (b)  the 
second  most  unstable  mode.  The  complex  frequencies  are  a>r/cur^ ax  —  0.77 
y/^max  **  0.09$  and  <t>r/comax  —  0.75,  y/c*>max  —  0.081  respectively.  The  growth  rates  are 
separated  more  than  for  the  previous  case  considered,  presumably  because  the  spatial  vari¬ 
ation  in  the  equilibrium  is  larger. 


FIG.  5.  Time  history  of  (a)  the  square  of  the  real  component  and  (b)  the  square  of  the 
total  perturbed  potential  for  a  fixed  location  xq. 
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FIG.  5b. 


FIG  6.  Mode  structures  obtained  from  simulations  showing  (a)  the  real  and  imaginary 
parts  and  (b)  the  square  of  the  mode  amplitude  as  functions  of  the  x  coordinate. 
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FIG.  6b. 
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FIG.  7.  Time  history  of  the  square  of  the  real  part  of  the  potential  at  a  fixed  location  x0 
as  a  function  of  time  for  the  case  vti/Ub  —  0.13. 


FIG.  8.  Test  panicle  time  history  in  the  kyAy  —  A x/Ln  plane  showing  the  spiral  orbits. 
Initially  the  particle  was  in  the  center  of  the  simulation  region.  This  is  from  the  simulation 
with  vu/Ub  —  0  13.  The  test  particle  initially  had  a  parallel  velocity  slightly  larger  than  the 
parallel  phase  velocity  of  the  wave. 
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FIG.  9.  Spatial  structure  of  the  potential  at  saturation  for  the  simulation  of  FIG  7  -  8 
The  contours  are  equally  spaced  with  the  dotted  lines  indicating  a  negative  value  and  the 
solid  line  indicating  a  positive  value. 
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FIG  10.  Test  particle  time  history  in  the  ky\y  -  A x/Ln  plane  showing  the  spiral  orbits. 
Initially  the  particle  was  in  the  center  of  the  simulation  region.  This  is  from  the  simulation 
with  v,,/  Ub  —  0.25.  The  test  particle  initially  had  a  parallel  velocity  slightly  larger  than  the 
parallel  phase  velocity  of  the  wave. 


FIG.  11.  Test  particle  motion  in  the  x  —  v.  plane  at  saturation  of  the  instability  for  the 
case  v,,/ Ub  —  0.19.  Here  all  of  the  particles  had  initial  x  locations  in  the  center  of  the  sys¬ 
tem.  Their  periodic  dimension  was  loaded  uniformly. 
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FIG  12.  Density  as  a  function  of  the  x  coordinate  at  saturation  of  the  field  energy  for  (a) 
the  background  ions  ,  (b)  the  stable  beam  ions  and  (c)  the  unstable  beam  ions.  Only  the 
unstable  beam  suffered  nonlinear  behavior.  The  parameters  of  Fig.  11  are  used  except 
v„/Ub  -  0.0. 
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FIG.  12b. 


FIG.  12c. 
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FIG.  13.  Density  as  a  function  of  the  x  coordinate  at  t  —  0  (  the  cosine  profile  )  and  at 
saturation  of  the  field  energy  for  (a)  the  background  ions  ,  (b)  the  stable  beam  ions  and 
(c)  the  unstable  beams  ions.  Both  the  unstable  beam  and  the  background  ions  suffer  non¬ 
linear  behavior.  The  saturation  value  was  e8<t>/Te  —  0.28. 
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FIG.  13c. 


FIG.  14.  Distribution  function  of  the  unstable  ion  beam  averaged  over  1/10  of  the  system 
length  at  t  —  0  ( the  smooth  curve  )  and  at  saturation  (a)  at  the  location  xj Lx  =  0.3  and 
(b)  at  the  location  xJLx  —  0.7. 
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FIG. 14b. 


Ret^.k.M 


FIG  15  Time  history  of  the  square  of  the  real  part  of  the  potential  at  a  fixed  location  x0. 
The  parameters  are  the  same  as  for  Fig.  13  except  that  the  0'h  and  1  were  kept  in  the 
periodic  direction  and  the  initial  potential  had  the  profile 
4>0(x)  -  (Te/e)(  0.5  +  Q.5cos(ttx/Lx)). 
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FIG.  16.  Density  as  a  function  of  the  x  coordinate  at  t  -  0  (  the  cosine  profile  )  and  at 
saturation  of  the  field  energy  for  (a)  the  background  ion  ions,  (b)  the  stable  ion  beam  and 
(c)  the  unstable  ion  beam.  The  parameters  are  the  same  as  for  Fig.  15.  The  saturation 
value  is  e8<t>/Te  —  0.13. 
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FIG.  16c. 


